home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / polynm.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  10.6 KB  |  301 lines  |  [TEXT/????]
  1. ;;; $Header: polynm.scm,v 1.9 88/09/01 22:01:31 GMT cph Exp $
  2. ;;; last modified 8/20/88 (mh)
  3.  
  4. ;;;;        Dense Univariate Polynomial Arithmetic
  5. ;;; Procedures for manipulating polynomials. A polynomial is represented as
  6. ;;;  a list of possibly-complex coefficients in ascending order, i.e.,
  7. ;;;    a0 + a1s + a2s^2 + . . . + ans^n ==> (a0 a1 a2 . . . an)
  8.  
  9. (if-mit
  10.  (declare (usual-integrations = + - * /
  11.                  zero? 1+ -1+
  12.                  ;; truncate round floor ceiling
  13.                  sqrt exp log sin cos)))
  14.  
  15. (define (add-2-polys p1 p2)
  16.   (cond ((null? p1) p2)
  17.     ((null? p2) p1)
  18.     (else (cons (+ (car p1) (car p2))
  19.                 (add-2-polys (cdr p1) (cdr p2))))))
  20.  
  21. (define add-polys (accumulation add-2-polys '()))
  22.  
  23. (define (sub-2-polys p1 p2)
  24.   (cond ((null? p1) (scale-poly -1 p2))
  25.     ((null? p2) p1)
  26.     (else (cons (- (car p1) (car p2))
  27.             (sub-2-polys (cdr p1) (cdr p2))))))
  28.  
  29. (define (sub-polys . polys)
  30.   (cond ((null? polys) '())
  31.     ((null? (cdr polys))
  32.      (scale-poly -1 (car polys)))
  33.     (else (sub-2-polys (car polys)
  34.                (apply add-polys (cdr polys))))))
  35.  
  36. (define (scale-poly factor poly)
  37.   (map (lambda (x) (* x factor)) poly))
  38.  
  39. (define (mul-2-polys p1 p2)
  40.   (if (or (null? p1) (null? p2))
  41.       '()
  42.       (add-2-polys (scale-poly (car p1) p2)
  43.          (cons 0 (mul-2-polys (cdr p1) p2)))))
  44.  
  45. (define mul-polys (accumulation mul-2-polys '(1)))
  46.  
  47. ;;; The following procedure returns a list of two lists, representing
  48. ;;; the quotient and the remainder polynomials.  For example:
  49. ;;;    (div-polys '(5 4 3 1) '(1 1)) ==> ((2 2 1) (3))
  50. ;;; corresponds to the division:
  51. ;;;
  52. ;;;     s^3 + 3s^2 + 4s + 5                      3
  53. ;;;    --------------------- = s^2 + 2s + 2 + -------
  54. ;;;           s + 1                            s + 1
  55.  
  56. (define (div-2-polys p1 p2)
  57.   (if (< (length p1) (length p2))
  58.       (list '() p1)
  59.       (let ((nr (reverse p1)) (dr (reverse p2)))
  60.     (let ((sf (/ (car nr) (car dr))))
  61.       (let ((dp (div-2-polys (reverse
  62.                   (sub-2-polys (cdr nr)
  63.                            (scale-poly sf (cdr dr))))
  64.                    p2)))
  65.         (cons (reverse (cons sf (reverse (car dp))))
  66.           (cdr dp)))))))
  67.  
  68. (define (div-polys . polys)
  69.   (cond ((null? polys) '(() (0)))
  70.     ((null? (cdr polys))
  71.      (list '() '(1)))
  72.     (else (div-2-polys (car polys)
  73.                (apply mul-polys (cdr polys))))))
  74.  
  75. (define (deriv-poly p) ; (a0 a1 ... aN) => (a1 2a2 ... NaN), (a0) => (0)
  76.   (let ((n (length p)))
  77.     (if (= n 1)
  78.         (list 0)
  79.         (generate-list (- n 1)
  80.                        (lambda (i)
  81.                          (let ((j (+ i 1)))
  82.                            (* j (list-ref p j))))))))
  83.  
  84.  
  85. (define (poly-value poly point)
  86.   (if (null? poly) ;Horner's method
  87.       0
  88.       (+ (car poly) (* point (poly-value (cdr poly) point)))))
  89.  
  90. ;;; The following procedure returns a list of the complex roots of a
  91. ;;; polynomial described by a list of its coefficients in the syntax
  92. ;;; (b0 b1 b2 . . . bn) = b0 + b1*s + b2*s^2 + . . . + bn*s^n.
  93. ;;; Coefficients may be arbitrary complex or real numbers.
  94. ;;; Initial-guess is an optional complex number argument; if
  95. ;;; omitted, the default value is '(.1 . 1).
  96. ;;; Example:
  97. ;;; ==>(poly->roots '(.0626457 .407983 .548892 1.41505 .574432 1))
  98. ;;; ((-5.48521e-2 . .965945)(-.14361 . .596993)(-.177508 . 0)
  99. ;;; (-.14361 . -.596993)(-5.48521e-2 . -.965945)) in about 15 sec.
  100. ;;; The control structure of this algorithm is essentially the same
  101. ;;; as that described by Cuthbert, p40, which in turn derives from
  102. ;;; a calculator program developed by H-P.
  103.  
  104. (define (poly->roots poly . initial-guess)
  105.   (define (find-a-root poly)
  106.     (define (poly-deriv poly point)
  107.       (let iter ((lst (cdr poly)) (k 1))
  108.     (if (null? lst)
  109.         0
  110.         (+ (* k (car lst))
  111.            (* point (iter (cdr lst) (1+ k)))))))
  112.  
  113.     (define (iter old-value-mag old-point new-point n m)
  114.       (let ((new-value (poly-value poly new-point)))
  115.     (let ((new-value-mag (magnitude new-value)))
  116.       (cond ((>= new-value-mag old-value-mag)
  117.          (cond ((< m 10)
  118.             (iter old-value-mag
  119.                   old-point
  120.                   (+ old-point (/ (- new-point old-point) 4))
  121.                   n
  122.                   (1+ m)))
  123.                ((< new-value-mag 1e-4) new-point)
  124.                (else (list 'step-size-too-small poly))))
  125.         ((< (magnitude (- new-point old-point)) 1e-5) new-point)
  126.         ((< n 50)
  127.          (iter new-value-mag
  128.                new-point
  129.                (- new-point
  130.                (/ new-value (poly-deriv poly new-point)))
  131.                (1+ n)
  132.                0))
  133.         (else (list 'no-root-found poly))))))
  134.     (let ((start-poly (poly-value poly initial-guess)))
  135.       (iter (magnitude start-poly)
  136.         initial-guess
  137.         (- initial-guess
  138.            (/ start-poly (poly-deriv poly initial-guess)))
  139.         0
  140.         0)))
  141.  
  142.   ;; continued on next page.
  143.  
  144.   (if (null? initial-guess)
  145.       (set! initial-guess (make-rectangular .1 1))
  146.       (set! initial-guess (car initial-guess)))
  147.   (if (< (length poly) 2)
  148.       '()
  149.       (let ((root (find-a-root poly)))
  150.     (if (number? root)
  151.         (cons (heuristic-round-complex root)
  152.           (poly->roots (car (div-polys poly (list (- root) 1)))
  153.                    initial-guess))
  154.         root))))
  155.  
  156. (define (roots->poly root-list)
  157.   (if (null? root-list)
  158.       (list 1)
  159.       (let ((poly-of-rest (roots->poly (cdr root-list))))
  160.     (sub-polys (cons 0 poly-of-rest)
  161.                (scale-poly (car root-list) poly-of-rest)))))
  162.  
  163. ;;; Create the zero-polynomial of unspecified degree
  164.  
  165. (define (the-zero-polynomial)
  166.   (list 0))
  167.  
  168. ;;; Given a polynomial p, return the procedure that computes p(x)
  169.  
  170. (define (poly->function p)
  171.   (lambda (x) (poly-value p x)))
  172.  
  173. ;;; Given a list of distinct abscissas xs = (x1 x2 ... xn) and a list
  174. ;;; of ordinates ys = (y1 y2 ... yn), return the Lagrange interpolation
  175. ;;; polynomial through the points (x1, y1), (x2, y2), ... (xn, yn).
  176.  
  177. (define (make-interp-poly ys xs)
  178.   ;; given a point list, return a poly that evalutates to 1 at the 
  179.   ;; first point and 0 at the others.
  180.   (define (unity-at-first-point point-list)
  181.     (let* ((x (car point-list))
  182.            (px (reduce * (map (lambda (u) (- x u)) 
  183.                      (cdr point-list)))))
  184.       (if (zero? px)
  185.           (error "MAKE-INTERP-POLY: abscissas not distinct"))
  186.           (scale-poly (/ 1 px) (roots->poly (cdr point-list)))))
  187.   (let loop ((p (the-zero-polynomial))
  188.              (points xs)
  189.              (values ys))
  190.     (if (null? values)
  191.         p
  192.         (let ((q (unity-at-first-point points)))
  193.           (loop (add-polys p (scale-poly (car values) q))
  194.                 (left-circular-shift points)
  195.                 (cdr values))))))
  196.  
  197. ;;; Given a function F, an interval [a, b], and a specified N > 0: we 
  198. ;;; generate a polynomial P of order N (and degree N-1) that interpolates 
  199. ;;; F at the "Chebyshev" points mapped onto [a, b]. We assume that the 
  200. ;;; absolute error function E(x) = |F(x) - P(x)| is unimodal between
  201. ;;; adjacent interpolation points. Thus E(x) has altogether N+1 "bumps"; 
  202. ;;; the largest of these has height Bmax and the smallest has height Bmin. 
  203. ;;; Of course Bmax is an error bound for the approximation of F by P on 
  204. ;;; [a, b]; but Bmin has heuristic significance as a lower-bound for the 
  205. ;;; reduced error that might be obtained by tuning the interpolation points 
  206. ;;; to meet the equiripple criterion. We return the list (P Bmax Bmin).
  207.  
  208. (define (get-poly-and-errors f a b n)
  209.   (let* ((c (/ (+ a b) 2))
  210.          (d (/ (- b a) 2))
  211.          (imap (lambda (x) (+ c (* d x)))) ;map [-1, 1] -> [a, b]
  212.          (points (map imap (cheb-root-list n)))
  213.          (p (make-interp-poly (map f points) points))
  214.          (abserr (lambda (x) (abs (- (f x) (poly-value p x)))))
  215.          (abserr-a (abserr a))
  216.          (abserr-b (abserr b))
  217.          (max-and-min-bumps
  218.            (let loop ((pts points)
  219.                       (bmax (max abserr-a abserr-b))
  220.                       (bmin (min abserr-a abserr-b)))
  221.              (if (< (length pts) 2)
  222.                  (list bmax bmin)
  223.                  (let ((x0 (car pts)) (x1 (cadr pts)))
  224.                    (let ((bump (cadr (brent-max abserr x0 x1 1e-6))))
  225.                      (cond ((> bump bmax) (loop (cdr pts) bump bmin))
  226.                            ((< bump bmin) (loop (cdr pts) bmax bump))
  227.                            (else (loop (cdr pts) bmax bmin)))))))))
  228.     (list p (car max-and-min-bumps) (cadr max-and-min-bumps))))
  229.  
  230. (define (get-poly-function-and-errors f a b n)
  231.   (let* ((c (/ (+ a b) 2))
  232.          (d (/ (- b a) 2))
  233.          (imap (lambda (x) (+ c (* d x)))) ;map [-1, 1] -> [a, b]
  234.          (points (map imap (cheb-root-list n)))
  235.          (p (lambda->procedure (lagrange (map f points) points)))
  236.          (abserr (lambda (x) (abs (- (f x) (p x)))))
  237.          (abserr-a (abserr a))
  238.          (abserr-b (abserr b))
  239.          (max-and-min-bumps
  240.            (let loop ((pts points)
  241.                       (bmax (max abserr-a abserr-b))
  242.                       (bmin (min abserr-a abserr-b)))
  243.              (if (< (length pts) 2)
  244.                  (list bmax bmin)
  245.                  (let ((x0 (car pts)) (x1 (cadr pts)))
  246.                    (let ((bump (cadr (brent-max abserr x0 x1 1e-6))))
  247.                      (cond ((> bump bmax) (loop (cdr pts) bump bmin))
  248.                            ((< bump bmin) (loop (cdr pts) bmax bump))
  249.                            (else (loop (cdr pts) bmax bmin)))))))))
  250.     (list p (car max-and-min-bumps) (cadr max-and-min-bumps))))
  251.  
  252. ;;; Given polynomial P(x), substitute x = r*y and compute the resulting
  253. ;;;  polynomial Q(y) = P(y*r).
  254.  
  255. (define (poly-arg-scale p r)
  256.   (if (zero? r)
  257.       (error "zero scale factor in POLY-ARG-SCALE"))
  258.   (let loop ((result '()) (factor 1) (poly p))
  259.     (if (null? poly)
  260.         (reverse result)
  261.         (loop (cons (* factor (car poly)) result)
  262.               (* factor r)
  263.               (cdr poly)))))
  264.  
  265.  
  266. ;;; Given polynomial P(x), substitute x = y+h and compute the resulting
  267. ;;;  polynomial Q(y) = P(y+h)
  268. ;;; Note: We use here the fact that DIV-polys returns '() for the quotient
  269. ;;;  if the degree of the denominator exceeds that of the numerator
  270.  
  271. (define (poly-arg-shift p h)
  272.   (let loop ((result '()) (poly p))
  273.     (if (null? poly)
  274.         (reverse result)
  275.         (let ((q (div-polys poly (list (- h) 1))))
  276.           (let ((quot (car q)) (rem (cadr q)))
  277.             (loop (cons (poly-value rem h) result)
  278.                   quot))))))
  279.  
  280.  
  281. ;;; Alter the coefficients of polynomial P so that its domain [a,b] is
  282. ;;;  mapped onto the canonical domain [-1,1].
  283.  
  284. (define (poly-domain->canonical p a b)
  285.   (if (<= b a)
  286.       (error "bad interval: must have a < b in POLY-DOMAIN->CANONICAL"))
  287.   (let ((c (/ (+ a b) 2)) (d (/ (- b a) 2)))
  288.     ;; p(x) [a,b] --> p(y+c) = q(y) [-d,d] --> q(d*z) = r(z) [-1,1]
  289.     (poly-arg-scale (poly-arg-shift p c) d)))
  290.  
  291.  
  292. ;;; Alter the coefficients of polynomial P so that its domain [-1,1]
  293. ;;;  is mapped onto [a,b].  This is the inverse operation to
  294. ;;;  POLY-DOMAIN->CANONICAL.
  295.  
  296. (define (poly-domain->general p a b)
  297.   (if (<= b a)
  298.       (error "bad interval: must have a < b in POLY-DOMAIN->GENERAL"))
  299.   (let ((c (/ (+ a b) 2)) (d (/ (- b a) 2)))
  300.     (poly-arg-shift (poly-arg-scale p (/ 1 d)) (- c))))
  301.